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We consider a plane periodical array of parallel cylindrical waveguides with evanescent coupling between 
, them. A new method for calculating of isofrequency curves based on the multiple Mie scattering formalism 
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(MMSF) is developed. This method is compared with the phenomenological model. The derivation of the 



^ ' phenomenological model by means of the MMSF is performed. The formulae for calculation of parameters 
•I— I 

^ of the phenomenological model are derived, such as propagation constants and coupling constant. 
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I. INTRODUCTION. 



Nowadays, much attention is devoted to periodical arrays of evanescently coupled optical waveg- 
uides. Such systems represent the particular case of low- dimensional photonic crystal structures. 
The general feature of such systems is the existence of photonic band structure that is analo- 
gous to the electron band structure in solids. Therefore some effects in photonic crystal structures 
may be analogous to some phenomena in solids jsl. 

In this paper we consider a plane array of parallel equidistant waveguides. We assume that the 
interaction between the waveguides is enough weak but not negligible. In this case, the eigenmodes 
of the array may be represented in a spirit of tight binding method taken from the solid state 
physics. It means that the eigenmodes of the array can be expressed in terms of the eigenmodes 
of the noninteracting waveguides. 

The eigenmodes for the j-th waveguide are described as follows j4| 

E,(r) = e-^-*+^''^-U,(a:„y,), 

(1) 

H,(r) = e— *+^^^^V,(x„y,), 

where w is a frequency of an eigenmode, the coordinates of a point r with respect 

to the axis of the waveguide, Pj is the propagation constant of the j-th waveguide. If Pj > u (the 
speed of light is assumed to be unit), the functions Uj(xj, \j{xj,yj) outside the waveguide 
decrease exponentially as the distance of the observation point from the waveguide increases. Thus, 
the mode is evanescent and it cannot be converted into a free photon. 

So, the eigenmodes of the array of weakly interacting waveguides may be represented in following 
way: 

N 



(2) 

N ^ ' 

If the distance between the waveguides is large enough, the coupling between only the nearest 
waveguides may be taken into account. Then, the equation for an eigenmode of the array reads 

'^^"^^ + ^'^'^""^ + 7(^,-i(^) + ^,+1 W) = 0, (3) 
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where 7 is a nearest neighbor couphng constant (for derivation see, for example, j4|). This 
equation or the analogous equations are usually used for simulation of optical effects in systems of 
interacting waveguides, such as optical Bloch oscillations js-?!, Zener tunneling js-lOl, dynamic 



localization 
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Il2|. Anderson localization jl3[. 
A principle drawback of Eq. is that the phenomenological constants Pj and 7 are unknown. 
They can be found from the experiment if one supposes that Eq. ([3]) is applicable to the system 
under study. 

However, in the important particular case of the cylindrical waveguides, one may propose the 
exact description of the optical properties of the array. In this case every eigenmode is characterized 
by the angular momentum m. The eigenmodes are described as follows 



H,„(r) = e-*-*+*"^<^^+*'^^™^ V,^(p,)- 



(4) 



Here m is the angular momentum, pj, (f)j are polar coordinates of a point r with respect to the 
axis of the waveguide, Ujm{pj), Vjm{pj) ~ Hm{>ijmPj), where Hm{x) is a Hankel function of the 



first kind, 



The rigorous formalism for description of array of cylindrical waveguides makes use of the exact 
solution for the electromagnetic wave scattering problem by an infinite cylinder. The proposed 
description is based on the possibility to generalize this solution for the case of many parallel 
cylinders — multiple Mie scattering formalism (MMSF) for the arrays of infinite cylinders 
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'his approach is similar to the multiple Mie scattering formalism for spherical particles 



22|. The MMS 



photonic crystals 
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can be used for investigation of scattering and transmission of light by 



17| or irregular systems of cylindrical waveguides |15|, for calculation of the 



eigenmode frequencies and band structures of a plane array of cylindrical waveguides [18] and of 
two-dimensional photonic crystals jl7|. 

In this paper we use the MMSF for calculation of the isofrequency curves for the array of 
identical waveguides. We consider the case when the waveguides of the array are situated close 
to each other and ascertain if the phenomenological approach is applicable for such a system. For 
this purpose we derive the phenomenological approach from the MMSF and develop the method 
for calculation the coupling constant. 



The MMSF is explained in Sect. [ITl In Sect. IIIII we discuss the connection between the MMSF 
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and the phenomenological approach and explain how the couphng constant 7 in Eq. ([3]) can be 
calculated. In Sect. IIVI we calculate the isofrequency curves of a plane array of infinite cylindri- 
cal waveguides. We compare the isofrequency curves calculated by means of MMSF with that 
calculated by the phenomenological model. In Conclusion we discuss the possibility of further 
development of method used in this paper. 

II. MULTIPLE MIE SCATTERING FORMALISM. 

Let us consider an array of parallel cylindrical dielectric rods. The axes of the rods are in 
the plane y = and they all are parallel to the 2;-axis. The refractive index of the j-th array is 
denoted as nj. Let the array is illuminated by the external field of the certain frequency uj and 
longitudinal wave vector /3: 

E^"* (r) = e-'^'+'^' E""* (x, y) , H'^"* (r) = 6"*'^*+*''" H""* (x, y) . (5) 
This field causes the response of the array. The field inside the j-th rod is 

E,(r) = Yl (^^- ^l^P^) - , 

(6) 

H, (r) = e-^-*+^^^ ^^'"'^^ (^^™ K^M + d,m Mi^^^(p,)) , < R. 

m 

Here pj, (pj, z are the cylindrical coordinates of the point r respectively to the axis of the j-th 
waveguide, and Uj = riju. The functions Mj^_^.^^(pj) and N^_,^^(pj) are linear superpositions of the 
Bessel functions. The partial amplitudes Cjm, djm determine the field inside the j-th rod. Below 
the factor e^"^*+*^^ is omitted, for short. 

The field scattered by the j-th rod may be represented in the form 

E,(r) = ^e^-*^ (a,^ M^^^(p,) - b,^ N^^„(p,)), 

(7) 

H,(r) = Y e^"^"^^ {^jm KM + ^rm Ml^^ip,)) , > R. 

m 

The functions M^^^(pj) and N^^^(pj) are linear superpositions of the Hankel functions of the 
first kind for the imaginary argument. 

On the other hand, the field Eq. ([7j) can be represented in the alternative form as an expansion 
in terms of functions M^^^(pi) and N^^^(p;) for any / 7^ j: 
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(8) 

H,(r) = e^^'^^ K^miPi) + Qrm ^UmiPi)) , ^ ^ J- 

m 

Let us emphasize that the Eqs. ([7]) and ([8]) represent the same field, i. e. the field scattered by 
the j-th waveguide. 



According to [I5l, one can relate the amplitudes p^j^, gj^ and aim, him as follows 

+ 0O +00 

n=— oo n=— oo 

where 

U'p^{u,(5) = K^.m{>ia-\3-l\)x\ ^ (10) 

[ (-1)™-" if /<j, 

X = a/cJ^"^-^ and Hm{x) is the Hankel function of the first kind. 
Let us introduce a notation 



E;.(r) = E EKr) = E E ^'"'^ (p;,nMi,„(p,) - g]^ Ni,^(p,)^ 
H;.(r) = 5^ HKr) = E E ^'"^^ (^U + im^l^miPj)) , p. > R- 



One can rewrite it in the form 



E;(r) = 5^ EKr) = E e'"""^' {p,m M-l^^ ' Qjm , 
H;(r) = 5^ HKr) = E e'""^' {pjm KpM + Qjm M^^m) , P. > ^• 



Here 



+ CX) 



l^j l^j n=-oo 

+ 00 

gjm = E^5- = E E U'j''Ju,P)bln. 
¥j I'T^j n=-oo 

Let us assume that the external field 



(11) 



(12) 



(13) 
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E-^r) = e'"""' {PL ^UM - QL K,m) , 

m 

H-^r) = e^'"^^ {PL KM + Q'm ^UnXP. 



(14) 



(15) 



Then, the field outside of the array may be represented in the form 

E(r) = E-*(r) + E,(r) + 5^EKr), 

H(r) = H-*(r) + H,(r) + 5^HKr), 
where the number j is arbitrary. 

The relations between fields outside and inside the j-th rod follow from the boundary conditions 
on its surface. These relations take the following form: 

=4".K/5) (16) 

"""'^ fm,{u3,P) (17) 

ymjj 

Taking into account Eq. ( IT3|) in Eq. (fT6|) one obtains the self-consistent system of equations 



N +00 



5-(u;, /3) r^i - E E /3) r i = • (18) 



1 ain^ 




rn 


[bin) 




[Qin 



bjm I l^j 71=— oo 

The system of equation Eq. flTSl) describes the response of the array on the external electromag- 
netic field, determined by the amplitudes P^, Qi^. At the same time, if one takes = Qi^ = 0, 
the Eq. (fTS!) describes the electromagnetic eigenmodes for the array under consideration. These 
modes are described as follows: 

\ A'' +00 



S-^i^.^) -E E U^^m{oo,f3) " =0. (19) 

\ bjm I l^j n=— oo y bin I 

The homogeneous system (fT9|) possesses a nontrivial solution only if 



det Sr^iu, /3) 6,i d^n - U^li^^, f3) = 0. (20) 
This equation allows to obtain the eigenvalues of /3 for the eigenmodes of the array. 
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In particular, for the single rod this equation takes the form 



det S-^iu,(3) =0. 



(21) 



The solutions of this equation fSj^ (depending on u) are the propagation constants of the j-th 
waveguide, that is assumed noninteracting with the other waveguides. One can see that these 
propagation constants are characterized by the angular momentum m, as it was mentioned above. 

Below we apply Eq. (fT8|) and Eq. (fT9|) to describe the optical properties of the array of the 
rods. 

III. RELATIONSHIP OF THE MULTIPLE SCATTERING FORMALISM AND THE 
PHENOMENOLOGICAL APPROACH. 

Let us derive the simplified equations which describe the optical properties of the array of the 
rods under consideration. 

Every rod is characterized by a set of its propagation constants /3jm, satisfying to Eq. f l2T]) . Let 
us notice that the propagation constants corresponding to the opposite angular momenta coincide. 



We suppose that the propagation constants of different waveguides differ slightly. Besides, the 



Therefore, we can consider the optical excitations originated from the propagation constants with 
fixed angular momentum m. Two cases are possible: m = and m ^ 0. For the first case 
one should take into account two partial amplitudes a^o, bjQ. For the second case the system of 
equations should include four partial amplitudes ajm, bjm and aj-m-, bj^-m, since the propagation 
constants for the angular momenta m and —m coincide. 

Below we take into account only the coupling between the nearest neighbors, since the coupling 
is evanescent. 

A. First case: m = 0. 

The first case is m = 0. In this case the main system of equations takes the form 



i. e. 




coupling Uj] 



{uj,f3) is weak and may be considered as a perturbation with respect to Sj^{u,(3). 
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bjo j i=j±i \biQ 



Below we omit the arguments cj, /9 for short. 
The matrix SJq is diagonal, 



S-o' =\^' ° I • (23) 



B, 



Here Aj, Bj are some functions of u and (3. So, Eq. (122|) separates in two independent systems 
of equations: 

A- «jo - 5^ U^o aio = 0, (24) 
i=j±i 

Bj bjo - ^io ^'0 = 0, (25) 
i=j±i 

Each propagation constant Pjq satisfy to one of the following equations: 

Aj{u,f3jo)=0, Bj{u,f3jo)=0. (26) 

Below we consider Eq. only, since for Eq. (12 5 p the derivation is similar. 
Since the coupling between the waveguides is weak, the isofrequency curve originating from any 
propagation constant is narrow. Therefore one can represent Aj in following way: 

A, = Df X - ^,o)- (27) 

Here f3jQ satisfies to the first of Eqs. (p6|) . 
Substituting this to Eq. (12^ . we obtain: 

jjlO 

(/3 - f3jo) a,o -Y. = 0. (28) 

i=j±i i 

Let us notice that U^q'^''^ = U^q^'^- Here we assume that ?7jg(w,/3) = t/j[](u;,/3jo) and that the 
value Uji^^'^{u, PjQ)/D^ is independent on j. So, introducing the notation 

j 



7 = (29) 



8 



we obtain 



(/3 - /3jo) ajo - 7 («i-i,o + «j+i,o) = 0. (30) 



The Eq. fl30|) possesses the nontrivial solutions only for eigenvalues of (3. 
The electric field outside the array of waveguides is 

N 

E{t, r) = * E ^''^ M^/3o(P.)- (31) 

The expression for the magnetic field is analogous. 

Here ^ means the sum over the eigenvalues of /3. We have added the argument /3 to the partial 
amplitudes ajo, bjQ, since the partial amplitudes depend on the eigenvalue p. 

Since the rods differ slightly and the interaction between them is weak, all the eigenvalues of /3 
are close to each other. So, one can suppose that M^^Q(pj) = M^^^,^jQ(pj). Let us introduce the 
notation 

%o(^) = ^'^^ «^o(/3). (32) 
So, the equation (13T|) takes the form: 

N 

E(t, r) = e--* ^Ip.M- (33) 

Taking into account Eq. ( !30|) . one can write the equation for ajo(-2): 



+ /^io j aio(^) + 7 («i-i,o(2) + «i+i,o(^)) = 0. (34) 



dz 

This equation coincides to Eq. ([3]). 

In the similar way one can derive the equation 

d 

where 



+ Pjoj bjo{z) + 7 (bj.i^oiz) + 6^+1,0(2:)) = 0. (35) 



7 = (36) 



and is determined by the equation 
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i?,=Df x(/3-/3,o). (37) 



B. Second case: m ^ 0. 



The second case is m 7^ 0. As it was mentioned above, we should take into account the partial 
amplitudes ajm, bjm and aj-m, bj-m both. So, the main system of equations takes the form 




/ aim 1 














(38) 



It is convenient to introduce the following notations: 



Tj3-^,rn _ jjj+l.m _ jjj-l-m _ jjj+l-m _ jj 
^ jm jm ^ j,~m ^ j,~m 

Tjj-l-m _ Tj-j+l-m _ Tjj-l,m _ Tjj+l,m _ y 
jm jn^ jy—'m, j,—m ^m' 



(39) 



Substituting this to Eq. ( l38l) . one gets 




0, 



(40) 



0. 



Below we show that there are two types of solutions. 

Let us suppose that the partial amplitudes aj^i bjm and aj_mi bj^^rn are connected by the 
following relation: 









\bj,-mj 




[bjml 



(41) 



Substituting this to (HOi) . one gets 
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(42) 



S-l^M hA - Y: (f/,n + V;„M-^) I =0. 

\bjm J l=j±l Y'ln 



Both equations in f H2|) should coincide. Therefore the matrix M should satisfy to following 
conditions: 

M-^ = M, 

(43) 

To find the possible forms of matrix M one should use a relation between the matrices Sjm and 
Sj^-^m- These matrices possess the form 



Sr^=\ . h Sr^-m=\ 1 : I , (44) 



iA C\ liA -C 

I ' ^j.—m I 

-C iBj \C iB 

where A, B and C are some real functions of u and /3. 

So, one can find easily, that there are only two possible forms of the matrix M: 

^ /i o\ . /-I o\ 

M = , or M = . (45) 

So, we see that the solutions of equation ( 138|) separates in two different types. 
For the solutions of the first type 

(^j,—m Ojjmi bj,—m ^jm? (46) 

and ajm, bjm satisfy the equation 

Sil - 5^ + ° ] (""A = 0. (47) 

\bjm I l=j±l \ Urn — V„i j \ / 

For the solutions of the second type 

^j,—m ^jmi bj,~m ^jm: (48) 
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and ajm, bjm satisfy the equation 



sil T''" - E r""'^"" ° If'"! =0. (49) 



bjm j l=j±\ 

Below we consider an equation 







1- 








I hrn j 





SJ^ r'-|=0, (50) 

bjm I l=j±l \blm 



where 



Um + Vm 

Urn '^m 



Wm = for the first case, 



(51) 



Wm = I I for the second case. 

Um + Vm 



Let Ujm = I ~ I be the solution of the equation 

bjm ^ 

Sjmi.^.Pjm) 'ajm = Q- (52) 

Remain that [3jm satisfies to the equation det S~^{uj, (3jm) = 0. The vector Uj^ is one of 
the two eigenvectors for matrix S~^{oj, f3jm) possessing a vanishing zero eigenvalue. Let Vjm be 
the other eigenfunction for the matrix Sjl^{u,f3jm), with fijm being the corresponding eigenvalue. 
Thus, 

jml'^' /^i™) ^jm — l^jm^jm- (^3) 

The vectors Ujm, ^jm are linearly independent. One assumes that uj^Ujm = vj^ Vj„j = 1. Let 
us find a solution of Eq. fl47p in the form 



^jm 
bjm 



-^jm U-jm ~l~ Bjm ^jm- (^^) 
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Since the coupling of the adjacent waveguides is a small perturbation to Sj^ in Eq. fl47p 



the vector | | is practically "parallel" to Ujm- For this reason, \Bjm\ ^ l^jml- Within the 

bjr 

perturbation approach, the value (3 — (3jm is a small parameter. Then, 

S-^{co, (3) ^ SJ^{UJ, fi^m) + (/3 - fiim) Djm, (55) 

where Djm is the derivative of the matrix Sj'^{u,/3) taken at the point /3 = f3jm- Substituting 
and to (130]) . one gets: 



— Wm (^AimUim + BimVlm^ =0. 

l=j±l 

With Eq. fl52|] being taken into account, the first order perturbation approach gives 



(56) 



f^jm BjjYi ^jm ~l~ (/3 Pjm) -^jm Dj^n ^jm ^ ^ -^Im ^jm 0. (57) 

l=j±l 

It is convenient to introduce a vector wjm completely defined by the conditions: 



^jm Ujm 1 j 



(58) 



Multiplying Eq. (157|) by wj^ results in the equation 



(/3 - f3j^) Ajra - Aira w]^ #„(/3) Uj>„ = 0, (59) 

l=j±l 

If the variation of f3jm is small as j changes, the variation of the product wj^ iym(/3) is 
small, as well. Therefore, one can neglect its dependence on j. Denoting 

7 = w]„iy„(/3)u,-„. (60) 

one obtains: 

{(5 - Pj^) Ajm - 1 + = 0. (61) 

The electric field outside the array is 
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N 



j=i /3 (62) 

+e-^'"'^^ (a,,_^(/3)M^^,_„(p,) -6,,_„(/3)N^^,_„(p,.))}- 
The expression for the magnetic field is analogous. Here the sum over /3 means the sum over 
the eigenvalues of longitudinal wave vector of the array. We have added the argument /3 to partial 
amplitudes since they may be different for different eigenmodes of the array. 
In the first approximation, 



Ajm{fi) Ujm, (63) 



where the upper sign is for the first case and the lower sign for the second case. 
The vector Uj^ doesn't depend on /3. The eigenvalues /3 differ slightly, so one can suppose that 
M^/3,±m(Pi) ~ M2^,„,±m(Pi)> ^l{i,±m{pj) ~ ^,„,±m(Pj ) • Introducing the notation 

one gets 

E(t,r) = e— * ^^™(^) {«^™ [f"^'' M^,^.^^(p,) ± M^,^,^_„(p,)) + 

i=i (65) 

+6,™ [f-^^ ^l,,^MTe—^^ N^,^_„(p,))}. 

From Eqs. (16T]) and it follows 

«^ + /3im^ ^jm(2) + 7 (^i-l,m(^) + Aj+i,^(z)) = 0. (66) 

This equation coincides to Eq. ([3]). 

IV. APPLICATION FOR ISOFREQUENCY CURVES CALCULATION. 

Consider an infinite array of identical waveguides. The optical eigenmodes in this system possess 
the form of Bloch waves: 
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E(t,r) = 6"^'^*+^^^+^^^ U(r), 

(67) 

H(t,r) = er'^'^+iP^+ik^ V(r), 

where U(r), V(r) are the periodical functions relatively to the coordinate x. Here k is the 
transverse quasi wave vector belonging to the interval — tt < k < tt (here the period of the array is 
assumed to be unit). For the fixed frequency u the longitudinal wave vector /3 is connected with 
the transverse quasi wave vector k, and the function P{k) is the so-called isofrequency curve. 

For the field outside the waveguides Eq. (E7I) results in the relations for the partial amphtudes 

a = a e'''^"- h- =h e^''^" ffiSl 

For the case of periodical array of identical waveguides the scattering matrices for all the 
waveguides are the same. Besides, the coupling coefficients [/j^(ci;,/3) depend on j — I. So the 
system fllSp takes the form 

^-^(a;,/3) hA + Y, E {I - j)a) hA =0. (69) 

\ bjm I l=^oo n \bln ) 



Substituting fl68l) to f l69l) . one obtains 



s;n\^,P) I ^"^ I + E ui{uj,p,k) 171=0, (70) 

where 




+ 00 

[/^(a;,/3,fc)= ^ Ul(uj,P, [l - j)a) e^'^C-^^ (71) 

l=—oo 

In the matrix form this system of equations can be written as 

L{uj,(3,k) x = 0, (72) 

where the matrix I/(/3, k) contains the scattering matrices Sl^^{u, (3) and the coupling coefficients 
U^{u, (3, k), and the column vector x contains the partial amplitudes am, bm- 
The nontrivial solutions of system (170|) exists when 



det L{u,/3,k) = 0. (73) 
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For fixed frequency u this equation implicitly determines the isofrequency curves P{k). 

The isofrequency curves can be derived also from the phenomenological model. Consider the 
equations (130|) and ( 16T]) . Assuming Ajm = e^^""^ Am and substituting this to (16T!) . one immediately 
gets an explicit expression for isofrequency curves: 

(3{k) = /3m + 2-f cos ka. (74) 

Here Pm is the propagation constant corresponding to the angular momenta m and — m (it is 
independent on the number j of a waveguide since all the waveguides are identical). The similar 
result can be obtained from Eq. (130|) after substitution ajo = e^^""^ Oq. 

As it was mentioned above, for m 7^ the system fl55]) possesses two types of solutions with two 
different coupling constants 7. So, the propagation constant 7^ gives rise to two different 

isofrequency curves. For the case m = the propagation constant gives rise to one isofrequency 
curve only. 

Here we compare the isofrequency curves calculated by means of the phenomenological model 
with the results of the rigorous model based on the multiple scattering formalism. I. e. we compare 
the results of calculations based on equations ( !73|) and ( fM|) . 

We take the array of waveguides of refractive index n^g = 1.554, and the refractive index of 
the medium outside the waveguides is nmed = 1.457. We suppose that the period of the array is 
unit, a = 1. The waveguides are supposed to be situated close to each other, i. e. the radii of the 
waveguides are R = 0.5. The velocity of light in vacuum is assumed c = 1. 

Below we use the multiple Mie scattering formalism for calculating several isofrequency curves 
originating from different propagation constants. The obtained isofrequency curves are compared 
with the prediction of the phenomenological model. To calculate the coupling constants 7 we use 
the formulae obtained in Sec. IIIII 

For angular momentum m = we take two propagation constants: /3i = 126.671 and /32 = 
126.704. The coupling constants for them are 71 = —3.92 x 10~^ and 72 = —3.78 x 10~^. 

The obtained isofrequency curves are presented in Fig. [T]and Fig. [2l The results of calculation 
by means of MMSF are presented by dots, and the predictions of the phenomenological model are 
shown by the solid curves. The horizontal lines show the propagation constants. (Here and below 
the isofrequency curves are plotted for < k < 7C since the function P{k) is even, /3(— fc) = P{k).) 
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One can see that the isofrequency curves obtained by MMSF and the phenomenological model 
almost coincide. 




Figure 1: Isofrequency curve originating from the propagation constant /? = 126.671 (m = 0). Dots for 
the curve obtained by MMSF, sohd Hue for the curve obtained by phenomenological model. 

For the angular momentum m = 1 we take two propagation constants also: fS^ = 131.099 
and /34 = 132.0092. For every of propagation constants two coupling constants exist. For /^s the 
coupling constants are 73 = 9.51 x 10~'^ and 73 = 5.97 x 10^^. For they are 74 = 7.90 x 10~^ 
and 7^' = 7.12 x 10"^ 

The obtained isofrequency curves are represented in Fig. [3] and Fig. |H One can see that the 
agreement between the results of MMSF and phenomenological model for the angular momentum 
m = 1 is much worth then for m = 0. In spite of this, the phenomenological model is applicable 
for the qualitative description of the isofrequency curves. 
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Figure 2: Isofrequency curve originating from the propagation constant /3 = 126.704 (m = 0). Dots for 
the curve obtained by MMSF, solid Hne for the curve obtained by phenomenological model. 




Figure 3: Isofrequency curves originating from the propagation constant P = 131.099 (m = 1). Dots for 
the curves obtained by MMSF, solid lines for the curves obtained by phenomenological model. 
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132,0075- 



132,0070 -I . 1 , 1 , 1 , 1 , 1 . 1 . 1 . 1 

7r/8 7i/4 371/8 nil Sti/S l>nlA Vn/S tt 

Transversal quasi wave vector k 



Figure 4: Isofrequency curves originating from the propagation constant (3 = 132.0092 (m = 1). Dots for 
the curves obtained by MMSF, solid lines for the curves obtained by phenomenological model. 

V. CONCLUSION. 

In this paper we considered the planar arrays of cylindrical rods by means of two methods: the 
phenomenological model and the multiple Mie scattering formalism. 

The MMSF has several advantages over the phenomenological method based on Eq. ([3]). First, 
the MMSF allows to calculate the behavior of the optical excitation for the case of the strong 
coupling between the waveguides, while the phenomenological method is applicable only for the 
case of the weak coupling. Second, the input data for the MMSF are the geometrical properties of 
the array and refractive indices of waveguides, while the phenomenological method requires some 
data that should be obtained experimentally, such as the propagation constants of waveguides and 
coupling constants. 

We demonstrated that for the case of evanescent coupling of rods the phenomenological model 
can be derived from MMSF. We developed the method to calculate the propagation constants 
[ijm and coupling constants 7. The applicability of the developed method is demonstrated for 
different isofrequency curves. The method represented in this work allows to produce the numerical 
simulation without need of experimental investigation of components of optical devices. 
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The method developed in this paper was used for isofrequency curves calculation for the case 
of weak interaction between the waveguides only. But it may be useful also for the systems with 
strong coupling between the waveguides. In this case the hybridization of modes with different 
angular momenta may take place due to the coupling. Mathematically it means that one can't 
neglect the coupling coefficients Uj^{u),/3) with n ^ m. In this situation the isofrequency curves 
may possess the shape much more complicated than the phenomenological model predicts. 

The MMSF represented in this paper is convenient only for the waveguides of cylindrical form, 
because in this case the scattering matrix can be calculated easily. However, this method can 
be applied for the waveguides of another shape, but in this case it would be more difficult 
to calculate the scattering matrix. Besides, the scattering by noncylindrical waveguides would 
mix the harmonics with different angular momenta. Mathematically it means that the scatter- 
ing matrix S{u!, (5) contains some "nondiagonal" elements describing the transition of harmonics 
e'^'t'i Mi^^(p,), e'^t>^ ^IpmiPj) to harmonics e'^t>^ M2^„(p,), e'^t>^ with n ^ m. Due to 

the existence of nonzero "nondiagonal" elements, the scattering matrix S{u!,/3) can't be separated 
to several matrices Sm{i^,P) with fixed m. 
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